M=1024;
num_batches=80;

% permeability of free space
u0= 4*pi*(10^-7);

save_dirc = '\';

%% Set up current sources

device_x_dimension = 2000*10^-6;    

%relative standard deviation of current
J_std=0.1;

%resolution
res=1024;

dx= device_x_dimension/res;

%padding
padding=20;

% number of pixels of image
number_x_pixels = res+2*padding;

im_x_dim=dx*number_x_pixels;

%standoff
standoff = 50*10^-6;
stadoff_std=10*10^-6;

%depth of current source
d=14*10^-6;

%current
tot_current=0.05;

%current density 
base_current=tot_current/(d*76e-6);

%optical diffraction/resolution
sigma=2*10^-6;

% x, y, and z of current source
x_vector = linspace(-im_x_dim/2,im_x_dim/2,number_x_pixels);
y_vector = linspace(-im_x_dim/2,im_x_dim/2,number_x_pixels);


% make meshgrids
[x_grid,y_grid] = meshgrid(x_vector,y_vector);
[XX,YY]=meshgrid((1:number_x_pixels)/number_x_pixels,(1:number_x_pixels)/number_x_pixels);

pixel_number_x = number_x_pixels;       
pixel_number_y = number_x_pixels;

list_kx = (2*pi/(im_x_dim))*((-pixel_number_x/2:pixel_number_x/2-1));
list_ky = (2*pi/(im_x_dim))*((-pixel_number_y/2:pixel_number_y/2-1));

% make sure no zeroes are in this list
kx = [];
ky = [];
for s=1:length(list_kx)
    if list_kx(s)~=0
        kx=[kx,list_kx(s)];
    elseif list_kx(s)==0
        kx=[kx,1e-9];
    end
end
for s=1:length(list_ky)
    if list_ky(s)~=0
        ky=[ky,list_ky(s)];
    elseif list_ky(s)==0
        ky=[ky,1e-9];
    end
end

% make matrices for kx and ky grid
[kx_matrix,ky_matrix] = meshgrid(kx,ky);

% calculate k perpendicular matrix
k_perp = sqrt(kx_matrix.^2+ky_matrix.^2);
    
XL=number_x_pixels;
YL=number_x_pixels;


z=zeros(1,M);
A=zeros(1,M);
Complex=zeros(1,M);
J_sign=zeros(1,M);
 

for nb=1:num_batches
    clear JDat_1024
    clear BDat_1024
    clear JDat_64
    clear BDat_64
    clear JDat_128
    clear BDat_128
    clear jx
    clear jy
    clear bx
    clear by
    clear bz
    clear BX
    clear BY
    clear BZ
    clear Dat_64
    clear Dat_128
    
    tic
    disp(nb)
    
    JDat_1024=zeros(M,1024,1024,2);
    BDat_1024=zeros(M,1024,1024,3);
    
    z=zeros(1,M);
    A=zeros(1,M);
    Complex=zeros(1,M);
    J_sign=zeros(1,M);
    
    
    %% Forward problem (generate magnetic field) Fourier calculations - slice by slice
    
    for k=1:M
        J=zeros(number_x_pixels,number_x_pixels);
        Jx=zeros(number_x_pixels,number_x_pixels);
        Jy=zeros(number_x_pixels,number_x_pixels);
        Bx=zeros(number_x_pixels,number_x_pixels);
        By=zeros(number_x_pixels,number_x_pixels);
        Bz=zeros(number_x_pixels,number_x_pixels);
        
        
        %Focusing on closer stand-off distance for Mark's Data
        z(k)=standoff+rand()*stadoff_std;
        
        %For each run, randomize whether current is positive or negative
        J_sign(k)=2*round(rand)-1;
        
        %Generating the current profile
        A(k)=randi(3);
        Complex(k)=randi(3);
        
        %Choose which calss of current you use
        if A(k)==1
            if Complex(k)== 1 
                [Jx_f,Jy_f]=RightAngle2(Jx,Jy,J_sign(k),XL,YL,base_current*random('Normal',1,0.33));
            elseif Complex(k) == 2
                [Jx_f1,Jy_f1]=RightAngle2(Jx,Jy,J_sign(k),XL,YL,base_current*random('Normal',1,0.33));
                [Jx_f2,Jy_f2]=RightAngle2(Jx,Jy,J_sign(k),XL,YL,base_current*random('Normal',1,0.33));
                Jx_f=Jx_f1+Jx_f2;
                Jy_f=Jy_f1+Jy_f2;
            elseif Complex(k) == 3
                [Jx_f1,Jy_f1]=RightAngle2(Jx,Jy,J_sign(k),XL,YL,base_current*random('Normal',1,0.33));
                [Jx_f2,Jy_f2]=RightAngle2(Jx,Jy,J_sign(k),XL,YL,base_current*random('Normal',1,0.33));
                [Jx_f3,Jy_f3]=RightAngle2(Jx,Jy,J_sign(k),XL,YL,base_current*random('Normal',1,0.33));
                Jx_f=Jx_f1+Jx_f2+Jx_f3;
                Jy_f=Jy_f1+Jy_f2+Jy_f3;
            end
        elseif A(k)==2
             if Complex(k)==  1 
                [Jx_f,Jy_f]=Slope2(Jx,Jy,J_sign(k),XL,YL,base_current*random('Normal',1,0.33));
            elseif Complex(k) == 2
                [Jx_f1,Jy_f1]=Slope2(Jx,Jy,J_sign(k),XL,YL,base_current*random('Normal',1,0.33));
                [Jx_f2,Jy_f2]=Slope2(Jx,Jy,J_sign(k),XL,YL,base_current*random('Normal',1,0.33));
                Jx_f=Jx_f1+Jx_f2;
                Jy_f=Jy_f1+Jy_f2;
            elseif Complex(k) == 3
                [Jx_f1,Jy_f1]=Slope2(Jx,Jy,J_sign(k),XL,YL,base_current*random('Normal',1,0.33));
                [Jx_f2,Jy_f2]=Slope2(Jx,Jy,J_sign(k),XL,YL,base_current*random('Normal',1,0.33));
                [Jx_f3,Jy_f3]=Slope2(Jx,Jy,J_sign(k),XL,YL,base_current*random('Normal',1,0.33));
                Jx_f=Jx_f1+Jx_f2+Jx_f3;
                Jy_f=Jy_f1+Jy_f2+Jy_f3;
            end
        elseif A(k)==3
             if Complex(k)==  1
                [Jx_f,Jy_f]=WideSlope2(Jx,Jy,J_sign(k),XL,YL,base_current*random('Normal',1,0.33));
            elseif Complex(k) == 2
                [Jx_f1,Jy_f1]=WideSlope2(Jx,Jy,J_sign(k),XL,YL,base_current*random('Normal',1,0.33));
                [Jx_f2,Jy_f2]=WideSlope2(Jx,Jy,J_sign(k),XL,YL,base_current*random('Normal',1,0.33));
                Jx_f=Jx_f1+Jx_f2;
                Jy_f=Jy_f1+Jy_f2;
            elseif Complex(k) == 3
                [Jx_f1,Jy_f1]=WideSlope2(Jx,Jy,J_sign(k),XL,YL,base_current*random('Normal',1,0.33));
                [Jx_f2,Jy_f2]=WideSlope2(Jx,Jy,J_sign(k),XL,YL,base_current*random('Normal',1,0.33));
                [Jx_f3,Jy_f3]=WideSlope2(Jx,Jy,J_sign(k),XL,YL,base_current*random('Normal',1,0.33));
                Jx_f=Jx_f1+Jx_f2+Jx_f3;
                Jy_f=Jy_f1+Jy_f2+Jy_f3;
            end
        end
        
        %random reflections and rotations
        
if randi(2)==2
    Jx_f=-Jx_f;
    Jy_f=-Jy_f;
end 

if randi(2)==2
    Jx_f=-rot90(rot90(Jx_f));
    Jy_f=-rot90(rot90(Jy_f));
end 

if randi(2)==2
    temp=Jx_f;
    Jx_f=-rot90(Jy_f);
    Jy_f=rot90(temp);
end 

        %Take 2D FFT of currents
        jx=fftshift(fft2(Jx_f));
        jy=fftshift(fft2(Jy_f));
        
        exp_factor = exp(-k_perp*(z(k)));
        
        bx = d*u0/2*(exp_factor.*(jy));
        by = d*u0/2*(exp_factor.*(-jx));
        %bz = d*u0/2*1i*((exp_factor.*(-ky_2D.*jx+kx_2D.*jy)./k_perp_2D));
        bz = d*u0/2*1i*((exp_factor.*(-ky_matrix.*jx+kx_matrix.*jy)./k_perp));

    
        %optical diffraction
        
        sf=sigma/2;
    
        PSF=exp(-.25*sf^2*(kx_matrix.^2+ky_matrix.^2));
    
        bx_blur=bx.*PSF;
        by_blur=by.*PSF;
        bz_blur=bz.*PSF;
        
        Bx(:,:) = real(ifft2(ifftshift(bx_blur)));
        By(:,:) = real(ifft2(ifftshift(by_blur)));
        Bz(:,:) = real(ifft2(ifftshift(bz_blur)));
        
        JDat_64(k,:,:,1)=BinImage(Jx_f(padding:end-padding,padding:end-padding),16);
        JDat_64(k,:,:,2)=BinImage(Jy_f(padding:end-padding,padding:end-padding),16);
        BDat_64(k,:,:,1)=BinImage(Bx(padding:end-padding,padding:end-padding),16);
        BDat_64(k,:,:,2)=BinImage(By(padding:end-padding,padding:end-padding),16);
        BDat_64(k,:,:,3)=BinImage(Bz(padding:end-padding,padding:end-padding),16);
        
%         JDat_128(k,:,:,1)=BinImage(Jx_f(padding:end-padding,padding:end-padding),8);
%         JDat_128(k,:,:,2)=BinImage(Jy_f(padding:end-padding,padding:end-padding),8);
%         BDat_128(k,:,:,1)=BinImage(Bx(padding:end-padding,padding:end-padding),8);
%         BDat_128(k,:,:,2)=BinImage(By(padding:end-padding,padding:end-padding),8);
%         BDat_128(k,:,:,3)=BinImage(Bz(padding:end-padding,padding:end-padding),8);
%         
        JDat_256(k,:,:,1)=BinImage(Jx_f(padding:end-padding,padding:end-padding),4);
        JDat_256(k,:,:,2)=BinImage(Jy_f(padding:end-padding,padding:end-padding),4);
        BDat_256(k,:,:,1)=BinImage(Bx(padding:end-padding,padding:end-padding),4);
        BDat_256(k,:,:,2)=BinImage(By(padding:end-padding,padding:end-padding),4);
        BDat_256(k,:,:,3)=BinImage(Bz(padding:end-padding,padding:end-padding),4);
        
%         JDat_512(k,:,:,1)=BinImage(Jx_f(padding:end-padding,padding:end-padding),2);
%         JDat_512(k,:,:,2)=BinImage(Jy_f(padding:end-padding,padding:end-padding),2);
%         BDat_512(k,:,:,1)=BinImage(Bx(padding:end-padding,padding:end-padding),2);
%         BDat_512(k,:,:,2)=BinImage(By(padding:end-padding,padding:end-padding),2);
%         BDat_512(k,:,:,3)=BinImage(Bz(padding:end-padding,padding:end-padding),2);

    toc
        
    end

    Dat_out_64.Jx=squeeze(JDat_64(:,:,:,1));
    Dat_out_64.Jy=squeeze(JDat_64(:,:,:,2));
    Dat_in_64.Bx=squeeze(BDat_64(:,:,:,1));
    Dat_in_64.By=squeeze(BDat_64(:,:,:,2));
    Dat_in_64.Bz=squeeze(BDat_64(:,:,:,3));
    Dat_in_64.z=z;
    Dat_in_64.Class=A;
    Dat_in_64.Numwire=Complex;
    
%     Dat_out_128.Jx=squeeze(JDat_128(:,:,:,1));
%     Dat_out_128.Jy=squeeze(JDat_128(:,:,:,2));
%     Dat__in_128.Bx=squeeze(BDat_128(:,:,:,1));
%     Dat_in_128.By=squeeze(BDat_128(:,:,:,2));
%     Dat_in_128.Bz=squeeze(BDat_128(:,:,:,3));
%     Dat_in_128.z=z;
%     Dat_in_128.Class=A; 
%     Dat_in_128.Numwire=Complex;
    
    Dat_out_256.Jx=squeeze(JDat_256(:,:,:,1));
    Dat_out_256.Jy=squeeze(JDat_256(:,:,:,2));
    Dat_in_256.Bx=squeeze(BDat_256(:,:,:,1));
    Dat_in_256.By=squeeze(BDat_256(:,:,:,2));
    Dat_in_256.Bz=squeeze(BDat_256(:,:,:,3));
    Dat_in_256.z=z;
    Dat_in_256.Class=A; 
    Dat_in_256.Numwire=Complex;
    Dat_in_256.device_x_dim=device_x_dimension;
%     
%     Dat_out_512.Jx=squeeze(JDat_512(:,:,:,1));
%     Dat_out_512.Jy=squeeze(JDat_512(:,:,:,2));
%     Dat_in_512.Bx=squeeze(BDat_512(:,:,:,1));
%     Dat_in_512.By=squeeze(BDat_512(:,:,:,2));
%     Dat_in_512.Bz=squeeze(BDat_512(:,:,:,3));
%     Dat_in_512.z=z;
%     Dat_in_512.Class=A; 
%     Dat_in_512.Numwire=Complex;
    
    
    
    %% Saving the Data
%     filename_out=[save_dirc,'TrainingDat_128_out_',num2str(nb,'%03.f'),'.mat'];
%     filename_in=[save_dirc,'TrainingDat_128_in_',num2str(nb,'%03.f'),'.mat'];
%     save(filename_out,'Dat_out_128','-v7.3')
%     save(filename_in,'Dat_in_128','-v7.3')
%     
    filename_out64=[save_dirc,'TrainingDat_64_out_',num2str(nb,'%03.f'),'.mat'];
    save(filename_out64,'Dat_out_64','-v7.3')
    filename_in64=[save_dirc,'TrainingDat_64_in_',num2str(nb,'%03.f'),'.mat'];
    save(filename_in64,'Dat_in_64','-v7.3')
    
%     filename_out=[save_dirc,'TrainingDat_512_out_',num2str(nb,'%03.f'),'.mat'];
%     filename_in=[save_dirc,'TrainingDat_512_in_',num2str(nb,'%03.f'),'.mat'];
%     save(filename_out,'Dat_out_512','-v7.3')
%     save(filename_in,'Dat_in_512','-v7.3')
%     
    filename_out=[save_dirc,'TrainingDat_256_out_',num2str(nb,'%03.f'),'.mat'];
    filename_in=[save_dirc,'TrainingDat_256_in_',num2str(nb,'%03.f'),'.mat'];
    save(filename_out,'Dat_out_256','-v7.3')
    save(filename_in,'Dat_in_256','-v7.3')
%    
end

%%

clf
colormap default 
k=15;
subplot(2,3,1)
imagesc(squeeze(JDat_64(k,:,:,1)))
title('Jx')
axis xy equal tight 
colorbar

subplot(2,3,2)
imagesc(squeeze(JDat_64(k,:,:,2)))
title('Jy')
axis xy equal tight 
colorbar

subplot(2,3,4)
imagesc(squeeze(BDat_64(k,:,:,1)))
title('Bx')
axis xy equal tight 
colorbar

subplot(2,3,5)
imagesc(squeeze(BDat_64(k,:,:,2)))
title('By')
axis xy equal tight 
colorbar

subplot(2,3,6)
imagesc(squeeze(BDat_64(k,:,:,3)))
title('Bz')
axis xy equal tight 
colorbar

